#######################################################################################
# Comparison of different IRT models on UN General Assembly voting data from 1970
# Author: Guillermo Rosas
# This file follows "CodeUNabstentions.R", written by G. Rosas 
#######################################################################################

# Require packages
library (MCMCpack)
library (MASS)
library (FastImputation)
library (superdiag)
library (coda)
library (mcmcplots)
library (sm)
library (foreign)
library (rjags)
library (runjags)
library (multicore)
library (random)
library (car)

source ("CodeUNabstentions.R") 
# When sourcing, make sure JAGS/MCMCpack statements in "CodeUNabstentions.R" are commented out,
# or be willing to wait a few days

# We compare regular CJR IRT against our best competing principals model
# Results are kept in MixChains (IRT) and MixChains3 (competing)
deviance.IRT <- MixChains[,grep("deviance", colnames(MixChains), fixed=TRUE)]
deviance.com <- MixChains4[,grep("deviance", colnames(MixChains4), fixed=TRUE)]
deviance.voeten <- MixChains5[,grep("deviance", colnames(MixChains5), fixed=TRUE)]
rbind (
  round (mean (deviance.IRT)),
  round (mean (deviance.com)),
  round (mean (deviance.voeten))
)
# Note that the deviance statistic for IRT is much lower than for competing,
# but this is only because the outcome variables are coded differently, sometimes as binary, sometimes as ordered

# Plot alphas and betas against each other, to see if there are many differences
# alpha=item difficulty
# beta =item discrimination
alpha.IRT <- MixChains[,grep("alpha", colnames(MixChains), fixed=TRUE)]
alpha.com <- MixChains4[,grep("alpha", colnames(MixChains4), fixed=TRUE)]

beta.IRT <- MixChains[,grep("beta", colnames(MixChains), fixed=TRUE)]
beta.com <- MixChains4[,grep("beta", colnames(MixChains4), fixed=TRUE)]

par (mar=c(3,3,1,1))
plot (c(-10,10),c(-10,10),type="n")
segments (x0=apply (beta.IRT, 2, quantile, 0.1), x1=apply (beta.IRT, 2, quantile, 0.9), y0=apply (beta.com, 2, quantile, 0.5), y1=apply (beta.com, 2, quantile, 0.5))
segments (x0=apply (beta.IRT, 2, quantile, 0.5), x1=apply (beta.IRT, 2, quantile, 0.5), y0=apply (beta.com, 2, quantile, 0.1), y1=apply (beta.com, 2, quantile, 0.9))
points (xy.coords (apply (beta.IRT, 2, quantile, 0.5), apply (beta.com, 2, quantile, 0.5)), pch=19)
abline (a=0, b=1, lty=2)

par (mar=c(3,3,1,1))
plot (c(-10,10),c(-10,10),type="n")
segments (x0=apply (alpha.IRT, 2, quantile, 0.1), x1=apply (alpha.IRT, 2, quantile, 0.9), y0=apply (alpha.com, 2, quantile, 0.5), y1=apply (alpha.com, 2, quantile, 0.5))
segments (x0=apply (alpha.IRT, 2, quantile, 0.5), x1=apply (alpha.IRT, 2, quantile, 0.5), y0=apply (alpha.com, 2, quantile, 0.1), y1=apply (alpha.com, 2, quantile, 0.9))
points (xy.coords (apply (alpha.IRT, 2, quantile, 0.5), apply (alpha.com, 2, quantile, 0.5)), pch=19)
abline (a=0, b=1, lty=2)
# For the most part, credible intervals are shorter for the competing principals model

# A few auxiliary objects
theta.IRT <- Theta.irt[new.ord,2]
theta.com <- Theta.compete.4[,2]
Beta.IRT <- colMeans (beta.IRT)
Beta.com <- colMeans (beta.com)
Alpha.IRT <- colMeans (alpha.IRT)
Alpha.com <- colMeans (alpha.com)
abstention.rate <- rowSums (is.na(rc.exp[,3:70]), na.rm=T) / ncol (rc.exp[,3:70])

# Do abstention rates differ by "principal"?
summary (abstention.rate[rc.exp$US.client==1])
summary (abstention.rate[rc.exp$Soviet.client==1])
summary (abstention.rate[rc.exp$US.client==0 & rc.exp$Soviet.client==0])

# Theta plots
ord.id <- order (theta.IRT)
par (mar=c(3,3,1,1))
plot (c(-4,4),c(-4,4),type="n")
segments (y0=Theta.compete.4[ord.id,1], y1=Theta.compete.4[ord.id,3], x0=Theta.irt[new.ord,2][ord.id], x1=Theta.irt[new.ord,2][ord.id])
segments (y0=Theta.compete.4[ord.id,2], y1=Theta.compete.4[ord.id,2], x0=Theta.irt[new.ord,1][ord.id], x1=Theta.irt[new.ord,3][ord.id])
text (xy.coords(Theta.irt[new.ord,2][ord.id], Theta.compete.4[,2][ord.id]), label=tolower(rownames(new.RC))[ord.id], cex=0.6)
text (xy.coords(Theta.irt[new.ord,2][ord.id], Theta.compete.4[,2][ord.id]-0.5), label=round (abstention.rate[new.ord][ord.id], 2), cex=0.6, col="red")
abline (a=0, b=1, lty=2)

# A better graph for presentation
# countries with obvious differences:
rel.position <- ifelse (Theta.compete.4[ord.id,3] < Theta.irt[new.ord,1][ord.id], "more left", ifelse (Theta.irt[new.ord,3][ord.id] < Theta.compete.4[ord.id,1], "more right", "same"))

# The following code produces Figure 1 in the AJPS article
par (mar=c(2,1,0,1), cex.axis=1, las=1)
# pdf ("UNrelativeThetaPositions.pdf", h=13, w=9)
plot ( c(-4,4), c(0,(length(theta.IRT))), type="n", ylab="", xlab="", main="", axes=F)
axis(1)
segments (x0=Theta.irt[new.ord,1][ord.id], x1=Theta.irt[new.ord,3][ord.id], y0=1:length(theta.IRT), y1=1:length(theta.IRT), col="gray")
points ( xy.coords ( theta.IRT[ord.id], 1:length(theta.IRT)), pch=19, cex=0.5, col="gray")
segments (x0=Theta.compete.4[ord.id,1], x1=Theta.compete.4[ord.id,3], y0=1:length(theta.IRT)+0.4, y1=1:length(theta.IRT)+0.4)
points ( xy.coords ( theta.com[ord.id], 1:length(theta.IRT)+0.4), pch=19, cex=0.5)
# segments (x0=Theta.sig.31[ord.ideal,1], x1=Theta.sig.31[ord.ideal,3], y0=1:98, y1=1:98, col="green")
# points ( xy.coords ( Theta.sig.11[ord.ideal,2], 1:98), pch=19, col="orange")
move.to.left  <- c("Trinidad and Tobago","Peru","Bolivia","Japan","Iceland","Brazil")
move.to.right <- c("Romania")
text (  xy.coords ( Theta.compete.4[ord.id,1][rel.position=="more left"], c(1:length(theta.com))[rel.position=="more left"]+0.5), pos=2, label=move.to.left, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][rel.position=="more right"], c(1:length(theta.com))[rel.position=="more right"]-0.5), label=move.to.right, pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][9], c(1:length(theta.com))[9]-0.5), label="Cuba", pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][123], c(1:length(theta.com))[123]-0.5), label="US", pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][5], c(1:length(theta.com))[5]-0.5), label="USSR", pos=4, cex=0.7)
# text (  xy.coords ( theta.com[ord.id]+1, 1:length(theta.com)), label=tolower(rownames(new.RC))[ord.id], cex=0.6  )
legend ("topleft", bty="n", pch=19, col=c("black","gray"), legend=c("Abstentions assumed MNAR (competing principals)","Abstentions assumed MAR"))
# dev.off()
# Countries with significantly different ideal points seem to be:
# Brazil, Iceland, Japan, Bolivia, Peru, Trinidad and Tobago (all US clients)
# Notable but less marked differences also appear in:
# Lesotho, Honduras, Argentina, Venezuela, Jamaica, Philippines (Lesotho not a US client)
# Most of these countries seem to be more to the left, once we consider the competing principals mechanism


# Comparison with Voeten's modesl
# Identify awkward countries in Voeten's first model
Voeten.country.order.2[new.ord][ord.id][c(14,17,124:126)]
y.country.order.2[new.ord][ord.id][c(14,17,124:126)]
y.country.order[ch.ord.voeten.2][new.ord][ord.id][c(14,17,124:126)]
# y.country.order.2 yields the same order as Voeten.country.order.2
# and is the result of y.country.order.2 <- y.country.order[ch.ord.voeten.2]


# The following code produces Figure C in the online appendix
par (mar=c(2,1,0,1), cex.axis=1, las=1)
# pdf ("UNrelativeThetaPositionsVoeten.pdf", h=13, w=9)
plot ( c(-4,4), c(0,(length(theta.IRT))), type="n", ylab="", xlab="", main="", axes=F)
axis(1)
# points ( xy.coords ( theta.IRT[ord.id], 1:length(theta.IRT)), cex=0.5, col="red")
segments (x0=Theta.compete.4[ord.id,1], x1=Theta.compete.4[ord.id,3], y0=1:length(theta.IRT)+0.4, y1=1:length(theta.IRT)+0.4)
points ( xy.coords ( theta.com[ord.id], 1:length(theta.IRT)+0.4), pch=19, cex=0.5)
segments (x0=Theta.voeten[ch.ord,1][new.ord][ord.id], x1=Theta.voeten[ch.ord,3][new.ord][ord.id], y0=1:126, y1=1:126, col="grey")
points ( xy.coords ( Theta.voeten[ch.ord,2][new.ord][ord.id], 1:126), pch=19, cex=0.5, col="gray")
move.to.left  <- c("Trinidad and Tobago","Peru","Bolivia","Japan","Iceland","Brazil")
move.to.right <- c("Romania")
voeten.left   <- c("Great Britain","South Africa","Portugal")
voeten.right  <- c("Albania", "Yemen")
# text (  xy.coords ( Theta.compete.4[ord.id,1][rel.position=="more left"], c(1:length(theta.com))[rel.position=="more left"]+0.5), pos=2, label=move.to.left, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][rel.position=="more right"], c(1:length(theta.com))[rel.position=="more right"]-0.5), label=move.to.right, pos=4, cex=0.7)
text (  xy.coords ( Theta.voeten[ch.ord,1][new.ord][ord.id][c(124:126)]-0.5, c(1:length(theta.com))[c(124:126)]+0.4), pos=2, label=voeten.left, cex=0.7)
text (  xy.coords ( Theta.voeten[ch.ord,1][new.ord][ord.id][c(14,17)]+0.2, c(1:length(theta.com))[c(14,17)]+0.3), label=voeten.right, pos=4, cex=0.7)

text (  xy.coords ( Theta.compete.4[ord.id,3][9], c(1:length(theta.com))[9]-0.5), label="Cuba", pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][123], c(1:length(theta.com))[123]-0.5), label="US", pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][5], c(1:length(theta.com))[5]-0.5), label="USSR", pos=4, cex=0.7)
# text (  xy.coords ( theta.com[ord.id]+1, 1:length(theta.com)), label=tolower(rownames(new.RC))[ord.id], cex=0.6  )
legend ("topleft", bty="n", pch=19, col=c("black","gray"), legend=c("Abstentions assumed MNAR (competing principals)","Abstentions as ordered outcomes (Voeten)"), cex=0.8)
# legend ("topleft", bty="n", pch=c(19,1,19), col=c("black","red","gray"), legend=c("Abstentions assumed MNAR (competing principals)","Abstentions assumed MAR","Abstentions imply weaker Nay (Voeten)"))
# dev.off()


# Now we try the same, but with Voeten's model that codes abstentions, and abstentions only, as 2s
par (mar=c(2,1,0,1), cex.axis=1, las=1)
# pdf ("UNrelativeThetaPositions.pdf", h=13, w=9)
plot ( c(-4,4), c(0,(length(theta.IRT))), type="n", ylab="", xlab="", main="", axes=F)
axis(1)
points ( xy.coords ( theta.IRT[ord.id], 1:length(theta.IRT)), cex=0.5, col="red")
segments (x0=Theta.compete.4[ord.id,1], x1=Theta.compete.4[ord.id,3], y0=1:length(theta.IRT)+0.4, y1=1:length(theta.IRT)+0.4)
points ( xy.coords ( theta.com[ord.id], 1:length(theta.IRT)+0.4), pch=19, cex=0.5)
segments (x0=Theta.voeten.2[ch.ord.voeten.2,1][new.ord][ord.id], x1=Theta.voeten.2[ch.ord.voeten.2,3][new.ord][ord.id], y0=1:126, y1=1:126, col="grey")
points ( xy.coords ( Theta.voeten.2[ch.ord.voeten.2,2][new.ord][ord.id], 1:126), pch=19, cex=0.5, col="gray")
move.to.left  <- c("Trinidad and Tobago","Peru","Bolivia","Japan","Iceland","Brazil")
move.to.right <- c("Romania")
voeten.left   <- c("Great Britain","South Africa","Portugal")
voeten.right  <- c("Albania", "Yemen")
text (  xy.coords ( Theta.compete.4[ord.id,1][rel.position=="more left"], c(1:length(theta.com))[rel.position=="more left"]+0.5), pos=2, label=move.to.left, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][rel.position=="more right"], c(1:length(theta.com))[rel.position=="more right"]-0.5), label=move.to.right, pos=4, cex=0.7)
text (  xy.coords ( Theta.voeten[ch.ord,1][new.ord][ord.id][c(124:126)]-0.5, c(1:length(theta.com))[c(124:126)]+0.4), pos=2, label=voeten.left, cex=0.7)
text (  xy.coords ( Theta.voeten[ch.ord,1][new.ord][ord.id][c(14,17)]+0.2, c(1:length(theta.com))[c(14,17)]+0.3), label=voeten.right, pos=4, cex=0.7)

text (  xy.coords ( Theta.compete.4[ord.id,3][9], c(1:length(theta.com))[9]-0.5), label="Cuba", pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][123], c(1:length(theta.com))[123]-0.5), label="US", pos=4, cex=0.7)
text (  xy.coords ( Theta.compete.4[ord.id,3][5], c(1:length(theta.com))[5]-0.5), label="USSR", pos=4, cex=0.7)
# text (  xy.coords ( theta.com[ord.id]+1, 1:length(theta.com)), label=tolower(rownames(new.RC))[ord.id], cex=0.6  )
legend ("topleft", bty="n", pch=c(19,1,19), col=c("black","red","gray"), legend=c("Abstentions assumed MNAR (competing principals)","Abstentions assumed MAR","Abstentions imply weaker Nay (Voeten)"))
# dev.off()





par (las=0)
cor (cbind (Theta.compete.4[,2], Theta.compete.bis[,2], Theta.compete[,2], Theta.irt[new.ord,2], Theta.voeten[ch.ord,2][new.ord]))
par (mar=c(3,3,1,1))
plot (c(-4,4),c(-4,4),type="n")
abline (a=0, b=1, col="red", lty=3)
segments (x0=Theta.compete.4[,1], x1=Theta.compete.4[,3], y0=Theta.voeten[ch.ord,2][new.ord], y1=Theta.voeten[ch.ord,2][new.ord], col="gray")
segments (x0=Theta.compete.4[,2], x1=Theta.compete.4[,2], y0=Theta.voeten[ch.ord,1][new.ord], y1=Theta.voeten[ch.ord,3][new.ord], col="gray")
text (xy.coords(Theta.compete.4[,2], Theta.voeten[ch.ord,2][new.ord]), label=tolower(rownames(new.RC)), cex=0.6)
mtext (side=1, line=2, text="Competing principals")
mtext (side=2, line=2, text="Abstentions as weaker Nay")
legend ("topleft", bty="n", col="black", legend=c(paste ("Competing principals deviance: ", round (mean (MixChains4[,ncol(MixChains4)]),2) , sep=""), paste("Abstentions imply weaker Nay deviance: ", round (mean (MixChains5[,ncol(MixChains5)]),2), sep="")))





# comparison of patterns of abstention of US satellites
bolivia <- unlist (RC[rownames(RC)=="BOLIVIA",])
iceland <- unlist (RC[rownames(RC)=="ICELAND",])
brazil  <- unlist (RC[rownames(RC)=="BRAZIL",])
france  <- unlist (RC[rownames(RC)=="FRANCE",])
unstate <- unlist (RC[rownames(RC)=="UNITED.STATES",])

table (bolivia, unstate, useNA="always")
table (iceland, unstate, useNA="always")
table (brazil,  unstate, useNA="always")
table (france,  unstate, useNA="always")

# Agreement indices
sum (diag (table (bolivia, unstate))) / sum (table (bolivia, unstate))
sum (diag (table (iceland, unstate))) / sum (table (iceland, unstate))
sum (diag (table (brazil, unstate))) / sum (table (brazil, unstate))
sum (diag (table (france, unstate))) / sum (table (france, unstate))


# summarize posterior probability of probabilities of voting given agreement, and disagreement
# These statistics appear in Table 4 in the AJPS article
round (colMeans (MixChains4[,263:266]), 2)
for (i in 263:266){
	print (round (quantile (MixChains4[,i], prob=c(0.1, 0.9)), 2))
}

# Calculate linear predictor for all country/vote pairs, both in IRT and Competing principals
pred.IRT <- pred.com <- matrix (NA, ncol=length(Beta.IRT), nrow=length(theta.IRT))
for (i in 1:length(theta.IRT)){
  for (j in 1:length(Beta.IRT)){
    pred.IRT[i,j] <- Beta.IRT[j]*theta.IRT[i]-Alpha.IRT[j]
    pred.com[i,j] <- Beta.com[j]*theta.com[i]-Alpha.com[j]
  }
}

prob.IRT <- pnorm (pred.IRT)
prob.com <- pnorm (pred.com)

sum ( abs (new.RC - prob.IRT), na.rm=T)
sum ( abs (new.RC - prob.com), na.rm=T)
sqrt (sum ( (new.RC - prob.IRT)^2, na.rm=T))
sqrt (sum ( (new.RC - prob.com)^2, na.rm=T))
# As expected, the sum of the squared residuals (actual vote minus predicted probability of vote) is much larger for competing principals
# this is also true for absolute deviations
# This is because competing principals also tries to "guess" at abstentions

# Number of vote choices predicted differently by IRT and competing principals, along with abstention rate
cbind (rowSums (abs (new.RC - round (prob.IRT)), na.rm=T ), rowSums (abs (new.RC - round (prob.com)), na.rm=T ), abstention.rate)
summary (glm (rowSums (abs (new.RC - round (prob.IRT)), na.rm=T ) ~ abstention.rate, family="poisson"))
summary (glm (rowSums (abs (new.RC - round (prob.com)), na.rm=T ) ~ abstention.rate, family="poisson"))
# Abstention rate is not driving at all the differences in prediction

# Now look at prediction discrepancies between prob.IRT and prob.com, but only for missing votes
# Totals
sum (abs (round(prob.IRT) - round(prob.com)))
sum (abs (round(prob.IRT) - round(prob.com))[!is.na(new.RC)])
sum (abs (round(prob.IRT) - round(prob.com))[is.na(new.RC)])

# Correct predictions from both models
sum (ifelse (new.RC==round (prob.IRT), 1, 0), na.rm=T)
sum (ifelse (new.RC==round (prob.com), 1, 0), na.rm=T)
sum (!is.na(new.RC))

# Percentages
sum (abs (round(prob.IRT) - round(prob.com))) / (nrow(new.RC)*ncol(new.RC))
sum (abs (round(prob.IRT) - round(prob.com))[!is.na(new.RC)]) / sum (!is.na(new.RC))
sum (abs (round(prob.IRT) - round(prob.com))[is.na(new.RC)]) / sum (is.na(new.RC))  # Again, the major difference is in the missing votes

# Maybe next step needs to incorporate uncertainty
# We have 1000 scans, too many to be practical; look at only 100
subsamp <- sort (sample (1:1000, 100))

th.irt <- MixChains[subsamp,c(1:126)]
th.irt <- th.irt[,new.ord]
th.com <- MixChains4[subsamp,c(1:126)]
bt.irt <- beta.IRT[subsamp,]
bt.com <- beta.com[subsamp,]
al.irt <- alpha.IRT[subsamp,]
al.com <- alpha.com[subsamp,]

pred.full.IRT <- pred.full.com <- array (NA, dim=c(100, length(theta.IRT), length(Beta.IRT)))
for (x in 1:100){
  for (i in 1:length(theta.IRT)){
    for (j in 1:length(Beta.IRT)){
      pred.full.IRT[x,i,j] <- bt.irt[x,j]*th.irt[x,i]-al.irt[x,j]
      pred.full.com[x,i,j] <- bt.com[x,j]*th.com[x,i]-al.com[x,j]
    }
  }
}

prob.full.IRT <- pnorm (pred.full.IRT)
prob.full.com <- pnorm (pred.full.com)

abs.irt <- apply (prob.full.IRT, 1, function (x) {sum (abs (new.RC - x), na.rm=T)})
abs.com <- apply (prob.full.com, 1, function (x) {sum (abs (new.RC - x), na.rm=T)})
sqr.irt <- apply (prob.full.IRT, 1, function (x) {sqrt(sum((new.RC - x)^2, na.rm=T))})
sqr.com <- apply (prob.full.com, 1, function (x) {sqrt(sum((new.RC - x)^2, na.rm=T))})

# Things do not change much once we account for uncertainty; squared and absolute Bayesian residuals are still larger for competing principals models
mean (abs.irt); mean (abs.com)
mean (sqr.irt); mean (sqr.com)

# Last thing to do, calculate a measure of uncertainty of each parameter
colSDs <- function (data) { apply (data, 2, sd)}
sd.th.irt <- colSDs (MixChains[,c(1:126)])
sd.th.com <- colSDs (MixChains4[,c(1:126)])
sd.bt.irt <- colSDs (beta.IRT)
sd.bt.com <- colSDs (beta.com)
sd.al.irt <- colSDs (alpha.IRT)
sd.al.com <- colSDs (alpha.com)
# Look at the means of these quantities: the uncertainty intervals (here, the SD of the posterior distributions of the parameters) are smaller for competing principals models (.com), but not by a huge lot

# Inspect what's happening with some of these votes. In the model with full recursivity there is way more uncertainty about missing votes of USSR and US
cbind (colMeans (Lead1.partial), colMeans (Lead1.full)) # Look for definitions of Lead1.partial and Lead1.full in "CJR_on_UN.R"

